Bulk and Collective Properties of a Dilute Fermi Gas in the BCS-BEC Crossover 
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We investigate the zero-temperature properties of a dilute two-component Fermi gas with attrac- 
tive interspecies interaction in the BCS-BEC crossover. We build an efficient parametrization of 
the energy per particle based on Monte Carlo data and asymptotic behavior. This parametriza- 
tion provides, in turn, analytical expressions for several bulk properties of the system such as the 
chemical potential, the pressure and the sound velocity. In addition, by using a time-dependent 
density functional approach, we determine the collective modes of the Fermi gas under harmonic 
confinement. The calculated collective frequencies are compared to experimental data on confined 
vapors of 6 Li atoms and with other theoretical predictions. 
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I. INTRODUCTION 



II. BULK PROPERTIES 



At zero temperature, the bulk energy per particle £ of 
a dilute Fermi gas can be written as 



A hot topic in nowadays' many-body physics is the 
study of Fermi gases at ultra-low temperature. Indeed, 
current experiments with atomic vapors are able to op- 
erate in the regime of deep Fermi degeneracy 0, Q ■ The 
two-component Fermi gases of these experiments are di- 
lute because the effective range Rq of the interaction is 
much smaller than the mean interparticle distance, i.e. 
UfRq "C 1 where kp — (S^n) 1 ^ 3 is the Fermi wave 
vector and n is the gas number density. Even in this 
dilute regime the s-wave scattering length a can be very 
large: the interaction parameter k F a can be varied over a 
very wide range using the Feshbach resonance technique, 
which permits to vary the magnitude and the sign of a. 
The available experimental data on 6 Li atoms are concen- 
trated across the resonance, where a goes from large neg- 
ative to large positive values 0, Q and where a crossover 
from a Bardeen-Cooper-Schrieffer (BCS) superfluid to a 
Bose-Einstein condensate (BEC) of molecular pairs has 
been predicted |E IE ■ 



In this paper we propose a reliable analytical fitting 
formula for the energy per particle of a homogeneous two- 
component Fermi gas, by analyzingthe fixed-node Monte 
Carlo data of Astrakharchik et al. 8] . From this analyti- 
cal formula it is straightforward to calculate several bulk 
properties of the system by means of standard themo- 
dynamical relations. This fitting formula enables us to 
calculate also the collective modes of the Fermi gas under 
harmonic confinement by using the hydrodynamic theory 
in the local-density approximation (LDA) , including also 
a quantum pressure term. We compare our results with 
other theoretical predictions [3, Hfl Ull LL2I H3, LLl and, 
in particular, with the experimental frequencies of the 
collective breathing modes 



£ = -e F e(x) 
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where e F = h 2 k F /(2m) is the Fermi energy and e(x) is a 
yet unknown function of the interaction parameter x = 
kpa. In the weakly attractive regime (—1 <C x < 0) one 
expects a BCS Fermi gas of weakly bound Cooper pairs. 
As the superfluid gap correction is exponentially small, 
the function e(x) should follow the Fermi-gas expansion 
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e(x) = l + -x- 



4(11- 2 In (2)) 
21^2 : 



(2) 



In the weak BEC regime (0 < x <C 1) one expects 
a weakly repulsive Bose gas of dimers. Such Bose- 
condensed molecules of mass tum = 2m and density 
njvf = njl interact with a positive scattering length 
a>M = 0.6a, as predicted by Petrov et al. ^(|. In this 
regime, after subtraction of the molecular binding energy, 
the function e(x) should follow the Bose-gas expansion 
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In the so-called unitarity limit (x = ±oo) one expects 
that the energy per particle is proportional to that of a 
non-interacting Fermi gas [Tsl ]. The fixed- node diffusion 
Monte-Carlo calculation of Astrakharchik et al. @ finds 



e(x = ±oo ) = 0.42 ±0.01 



(4) 



while an analogous calculaton of Carlson et al. [jjj 
gave e(x = ±oo) = 0.44 ± 0.01. The calculation of As- 
trakharchik et al. H is quite complete and gives the be- 
havior of the energy of system across the unitarity limit. 
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BCS (y < 0) BEC (y > 0) 

expression value expression value 
ai £ 0.4200 £ 
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0.4200 
0.3692 
1.0440 
1.4328 



-CH2 - 



[fitted] 



0.4200 
0.2674 
5.0400 
0.1126 



fa a 2 Q 3 /3i/C- 0.5523 q 2 q 3 /3i/C+ 0.4552 



TABLE I: Parameters of the fitting function (6). 



It is a standard convention to use the inverse interac- 
tion parameter y = X/x = l/(kpa) as the independent 
variable. In Fig. 1 we plot the data of e(y) reported by 
Astrakharchik et al. 0, |2(J . On the basis of the data of 
Carlson et al. 0] , Bulgac and Bertsch 01 proposed the 
following behavior of e(y) near y = 0: 



e(y) = £-Cy + 



(5) 



with £ = 0.44 and ( = 1 for both positive and negative y. 
The denser data of Ref. [8j suggest instead a continuous 
but not differentiable behavior of e(y) near ?/ — 0, namely 
with £ = 0.42 and ( = £_ = 1 in the BCS region (y < 
0) but ( = (+ = 1/3 in the BEC region (y > 0). As 
expected, for large |y|, the Monte-Carlo data shown in 
Fig. 1 follow the asymptotic trends of Eq. (1) and Eq. (2). 

We propose here the following analytical fitting for- 
mula 



e(y) — ot\ — c*2 arctan az y 



Pi 



\y\ 



(6) 



interpolating the Monte Carlo energy per particle and 
the limiting behaviors for large and small \y\. Here the 
parameter ol\ is fixed by the value £ of e(y) at y = 0, 
the parameter 012 is fixed by the value of e(y) at y = 00, 
and 03 is fixed by the asymptotic 1 jy coefficient of e(y) at 
large \y\ (Eqs. (2) and (3)). The ratio P1/P2 is determined 
by the linear behavior £ of e(y) near y = 0. The value 
of 0i is then determined by minimizing the mean square 
deviation from the Monte-Carlo data j2j). Of course, we 
consider two different set of parameters: one set in the 
BCS region (y < 0) and a separate set in the BEC region 
(y > 0). Table 1 reports the values of these parameters. 

Figure 1 compares this fitting function (solid curve) 
to the Monte Carlo data. For the sake of completeness, 
in Fig. 1 we also show the dotted curve obtained with 
the [2,2] Pade approximation of Kim and Zubarev [Tlj . 
based only on the asymptotes and the Monte-Carlo value 
|19|] at y — 0. Our parametric formula is more accurate, 
especially around y = 0. 

The advantage of a functional parametrization of e(y) 
is that it allows straightforward analytical calculations of 
several ground-state physical properties of the bulk Fermi 
gas [2^. For example, the chemical potential fj, is given 
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FIG. 1: The energy per particle e(y) as defined in Eq. (1), 
where y = 1/x = l/(fejra). Solid circles represent the fixed- 
node Monte-Carlo data of Ref. 0- Solid line is the parametric 
function (6) based on the values of Table 1. Dotted line is the 
Pade approximation of Ref. [Tll| . Dashed lines represent the 
asymptotic expressions Eq. (2) and Eq. (3). 
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FIG. 2: Chemical potential fi/ep, pressure P/(nei?) and 
sound velocity c 3 /vf, obtained from the parametric func- 
tion (6). Dashed line represents the simple model m/ef = 
I - I arctan(7T2//2) of Ref. [lj. 
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by 
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dn 



e F (e(y) - |e'(y)) , (7) 



as found by using Eq. (1) and taking into account that 
dy/dn = —y/(3n), while the pressure P reads 



(8) 



The sound velocity c s is instead obtained as c 2 
(n/m)du/dn, from which we get 
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y r>(y) + y ^"(y) 
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where vp — (2cf /m) 1 / 2 is the Fermi velocity. Figure 2 
reports the chemical potential /i, the pressure P and the 
sound velocity c a as a function of y. Our theory predicts 
that all these macroscopic properties show a kink at the 
unitarity point, due to £- ^ C+- Figure 2 shows also the 
curve of the chemical potential obtained with the simple 
analytical model proposed by Combescot and Leyronas 
[l2^ . The sound velocity c s is accessible experimentally, 
and the dotted curve of Fig. 2 is our prediction of the 
way c s evolves from to zero through the BCS- 

BEC crossover. 



III. HARMONICALLY CONFINED GAS 

We consider now the effect of confinement due to an 
external anisotropic harmonic potential 



U(r) 



I 2 2 1 2 2\ 
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where uj p is the cylindric radial frequency and lo z is the 
cylindric longitudinal frequency. Assuming that the den- 
sity field n(r, t) varies sufficiently slowly (this assumption 
is at the basis of the LDA), at each point r the gas can 
be considered in local equilibrium, and the local chemical 
potential is n[n(r, t)]. Within the LDA, the dynamics can 
be described by means of the hydrodynamic equations of 
superfluids 



9n „ , . 
- + V.(nv) 



0, (11) 



On ( 1 \ 

m^+V U[n(r 1 t)] + U(r) + -mv 2 ) = 0, (12) 

where v(r, t) is the velocity field, and /j,[n] is the chem- 
ical potential of Eq. (7). It has been shown by Cozzini 
and Stringari [23| that assuming a power-law dependence 
H = hq n 1 for the chemical potential (polytropic equa- 
tion of state horn Eqs. (11-12) one finds analytic 
expressions for the collective frequencies. In particular, 
for the very elongated cigar-shaped traps used in recent 
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FIG. 3: Frequencies of the longitudinal and radial collective 
breathing modes of the 2-component Fermi gas as a function 
of a H /(N 1/6 a) = 2 1/2 3 1/6 /(k F a) = 2 1/2 3 1/6 y. Sl z experi- 
mental data from Ref. 0], and fi p experimental data from 
Ref. 0a (squares) and Ref. Qb (circles), both for 6 Li. Joined 
dots with solid lines: TDDFT of Eq. (16), with N = 4 x 10 5 , 
lu p /lo z = 31.9 for Q z and with N = 2 x 10 5 , uj p /uj z = 22.1 for 
Q p . Dot-dashed line: variational approximation to TDDFT, 
Ref. [Tll |. Dashed line: polytropic hydrodynamic approxima- 
tion Eq. (13-14). Dotted line: the mean-field BCS of Ref. 0. 



experiments {uj r /oj z > 20), the collective radial breathing 
mode frequency f2 p is given by p3 | 



a 



v/2( 7 +l)^ p , 



(13) 



while the collective longitudinal breathing mode Q z is 



? \l UJr, 

7 + 1 



(14) 



In our problem we introduce an effective polytropic 
index 7 as the logarithmic derivative of the chemical po- 
tential fx, that is 
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\e(y)- 2 fe'(y) + ^e"(y) 
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(15) 



We have verified that indeed 7 remains relatively close to 
unity for all y: the results of the local polytropic equation 
are thus useful to have a simple analytical prediction of 
the collective frequencies. Based on this polytropic hy- 
drodynamic approximation (PHA), by using Eq. (6) we 
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obtain the breathing-mode frequencies shown in Fig. 3 as 
dashed lines. 

The analytical prediction of Eqs. (13-15) can be im- 
proved by releasing the polytropic approximation and 
explicitly integrating Eqs. (11-12). We have done such 
a calculation by including also a quantum pressure term 
(— h 2 V 2 y/n) / (2myfn) in Eq. (12). In practice, one must 
solve the following time-dependent nonlinear Schrodinger 
equation 



ih—ip(T,t) = 



.—W 2 + U(r) + ^n(r,t)} 



where ip(r,t) is the superfluid wave function such that 
n(r,t) = |^>(r, t)| 2 , v = hVhi{ip/n 1 l' 1 )/{i'm) 1 and n[n] 
is the chemical potential of Eq. (7). Equation (16) 
can be interpreted as the Euler-Lagrange equation of a 
time-dependent density functional theory (TDDFT) . 
In Ref. 0|, Eq. (16) is approached via a variational 
scheme. Here instead, Eq. (16) is solved numerically 
by using a finite-difference Crank-Nicolson predictor- 
corrector scheme [24] . First we obtain the ground state 
by integrating Eq. (16) in imaginary time. Then we let 
a slightly perturbed wave function evolve in real time 
for approximately one period of oscillation of the low- 
est (longitudinal) frequency fl z . In the same time span, 
the density also undergoes several radial oscillations of 
frequency f2 p . We extract both frequencies by fitting the 
mean square widths of n(r, t) with the sum of two cosines. 
The breathing-mode frequencies obtained in this way are 
shown in Fig. 3 as dots joined by solid lines. 

The quantum pressure term is important for small 
number of atoms, as it improves the determination of 
the density profile close to the surface of the vapor cloud 
|ll| . For the number of particles of the experiments 
(N > 10 5 ), the quantum pressure term is a relatively 
small correction, and according to our calculation it ac- 
counts for about 0.5% of the total energy. 

Figure 3 shows substantial accord between PHA and 
TDDFT. The main difference is the location of the pre- 
dicted maximum in the bosonic region y > 0. The dif- 
ferences are not only due to the approximations involved 
in the PHA, but also to numerical errors in TDDFT in- 
troduced by space and time discretization (estimated to 
less than 1%)|25|. We observe that, according to our 
predictions, the collective frequencies reach the asymp- 
totic large- 1 y | limits more slowly than the theory of Hu et 
al. 0] based on mean-field BCS Bogoliubov-de Gennes 
equations within LDA. In particular, in the BEC region, 
our numerical and analytical results (see the inset of 
Fig. 3) show that, contrary to the mean-field prediction, 
rip approaches its asymptotic value 2 lo p passing through 
a local maximum. il z has the same nonmonotonic be- 
havior while reaching a/5/2 lo z for large y. This qualita- 
tive behavior was previoulsy suggested by Stringari 
The different asymptotic behavior of the BCS mean-field 
frequencies is due to the neglect of beyond mean-field 
corrections. Note also that further discrepancies on the 
BEC side are due to the mean-field relation a rn = 2a s 



rather than a m = 0.6a s as provided by four-body scat- 
tering and used in our calculation. The Q p curve 
computed in Ref. HU (dot-dashed line in Fig. 3) agrees 
rather well with our calculation. 

Different theories are compared to the experimental 
data by Kinast et al. ylfor the radial mode f2 p and those 
by Bartenstein et al. [4( for the longitudinal mode Q z . 
In Fig. 3, we use the standard variable aH/{N x ' 6 a) = 
2 1 / 2 3 1 / 6 y and follow Ref. (2|| to determine the scattering 
length a as a function of the magnetic field B near the 
Feshbach resonance: 



a b [1 + a(B - B )} 



B-B 



(17) 



where B a = 83.4149 mT, a b = -1405 a , A = 30.0 mT, 
and a = 0.0040 (mT) -1 . For the longitudinal frequency 
Cl z , our results are in quantitative agreement with the 
experimental data, not unlike the mean-field prediction. 
The accord is less good for the radial mode. Experimen- 
tal uncertainty of the position of the resonant field Ba 
could partly account for these discrepancies. In particu- 
lar the upward feature near a^/(A rl / 6 a) ~ —1.5 has been 
related 0, 0] to the breaking of the Cooper pairs (due 
to sizeable ratio between the collective energy hTt p and 
the gap energy causing a failure of the hydro- 

dynamical approximation. Finite-temperature and non- 
LDA effects not taken into account in the theories could 
be relevant. Note also that the experimental situation 
is not completely clear. The experimental measurement 
of f2 p performed by Bartenstein et al. £| (not shown in 
Fig. 3) disagree with the data of Kinast et al. @. In 
particular, Bartenstein et al. Q find VL p /uo p ~ 1.6 at 
the unitarity limit y = 0, instead of the expected value 
VL p /w p = v/lO/3 = 1.82 obtained from Eq. (13) for 7 = 2, 
characteristic of the (renormalized) free Fermi gas. 



IV. DISCUSSION 

We propose analytic expressions for the equations of 
state of a uniform dilute Fermi gas across the BCS-BEC 
transition. These expressions are based on recent Monte- 
Carlo data and well-established asymptotic expansions. 
By using a hydrodynamic local-density approximation we 
include the effect of harmonic confinement. We compare 
the predictions of this approach with experimental fre- 
quencies of confined 6 Li vapors. Other predicted physi- 
cal quantites can be accessed by future experiments. The 
hydrodynamic approach is improved, to address small 
number of atoms, by including a quantum-pressure term. 
Indeed, our parametric formula (6) provides an accurate 
expression for the fi term in Eq. (16) through Eq. (7), 
which gives a better determination of the density and 
of the collective frequencies. This generalized quantum 
hydrodynamic approach, which takes into account the 
beyond-mean-field corrections, provides a reliable tool to 
determine the density profile of the Fermionic cloud and 
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to investigate its collective dynamical properties, includ- to new experimental determinations of Q p . 
ing also mode coupling and anharmonic oscillations. 
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